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ABSTRACT 

Collapse of massive stars may result in formation of accreting black holes in their inte- 
rior. The accreting stellar matter may advect substantial magnetic flux onto the black 
hole and promote release of its rotational energy via magnetic stresses (the Blandford- 
Znajek mechanism). In this paper we explore whether this process can explain the 
stellar explosions and relativistic jets associated with long Gamma-ray-bursts. In par- 
ticularly, we show that the Blandford-Znajek mechanism is activated when the rest 
mass-energy density of matter drops below the energy density of magnetic field in the 
very vicinity of the black hole (within its ergosphere). We also discuss whether such 
a strong magnetic field is in conflict with the rapid rotation of stellar core required in 
the collapsar model and suggest that the conflict can be avoided if the progenitor star 
is a component of close binary. In this case the stellar rotation can be sustained via 
spin-orbital interaction. In an alternative scenario the magnetic field is generated in 
the accretion disk but in this case the magnetic flux through the black hole ergosphere 
is not expected to be sufficiently high to explain the energetics of hypernovae by the 
BZ mechanism alone. However, this energy deficit can be recovered via additional 
power provided by the disk. 

Key words: black hole physics - supernovae: general - gamma-rays: bursts - meth- 
ods: numerical - MHD - relativity 



1 INTRODUCTION 

The most popular model for the central engines of long 
Gamma-ray-burst (GRB) jets is based on the "failed su- 
pernova" scenario of stellar collapse, or "collapsar", where 
the iron core of progenitor star forms a BH (Wo osley|1993 1. 
If the specific angular momentum in the equatorial part of 
progenitor exceeds that of the marginally bound orbit of the 
BH then the collapse becomes highly anisotropic. While in 
the polar region it may proceed more or less uninhibited, at 
least for a while, the equatorial layers form dense and mas- 
sive accretion disk. The gravitational energy released in this 
disk can be very large, more then sufficient to overturn the 



collapse of outer layers and drive GRB outflows (MacFadyen 
& Wooslcy 19991. A similar configuration can be produced 
via inspiraling of a BH or a neutron star into the compan- 
ion star during the common envelope phase of a close binary 



( Tutukoy fc Iungelson|1979l|Frayer fc Woosley|1998| |Zhang 



& Fryer 20011. In any case, the observed duration of long 



GRBs, between 2 and 1000 seconds, imposes a strong con- 
straint on the size of progenitors - their light crossing time 
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should be significantly shorter then the burst duration. This 
implies that the progenitors must be compact stars stripped 
of their hydrogen envelopes. This conclusion agrees with the 
absence of hydrogen lines in the spectra of supernovae iden- 
tified with GRBs (Moreover, since some of these supernovae 
are of Type Ic their progenitors may have lost their helium 
envelopes as well.). Massive solitary stars can lose their en- 
velopes via strong stellar winds. However, this may lead to 
unacceptably low rotation rates by the time of stellar col- 
lapse ( Heger et al.||2005 |. In the case of binary system the 
loss of envelope can be caused by the gravitational interac- 
tion with companion. 

The currently most popular mechanism of powering 
GRB jets is the heating via annihilation of neutrinos and 



anti-neutrinos produced in the disk ( MacFadyen & Woosley 



1 1999p . The energy deposited in this way is a strong func- 
tion of the mass accretion rate which must be above ~ 
O.O5M0S" 1 in order to agree with the observational con- 
straints on the energetics of GRBs jP opham et al.||l999[ ). 
Such high accretion rates can be provided in the collapsar 
scenario and the numerical simulations by MacFad yen fc| 
|Woosley| ( |1999[ ) and |Aloy et al."| (pOOO!) have demonstrated 
that sufficiently large energy deposition in the polar region 
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above the disk may indeed result in fast collimated jets. 
However, we have to wait for simulations with proper im- 
plementation of neutrino transport before making final con- 
clusion on the suitability of this model - the long and com- 
plicated history of numerical studies of neutrino-driven su- 
pernova explosions teaches us to be cautious. The first at- 
tempt to include neutrino transport in collapsar simulations 
was made by Nagataki et al. ( 2007[) and they did not see 



neutrino-driven polar jets - the heating due to neutrino an- 
nihilation could not overcome the cooling due to neutrino 
emission]^] In any case, it is unlikely that this model can ex- 
plain the bursts that are longer than 100s as by this time 
the mass accretion rate is expected to drop significantly. 
For example the accretion rate in the simulations by |Mac-1 
Fadyen & Woosley ( 1999 1 has already drifted below the crit- 
ical M ~ O.O5M s"Mj}T t = 20s. 

The two alternatives to the neutrino mechanism which 
are also frequently mentioned in connection with GRBs are 
the magnetic braking of the accretion disk (|Blandford & 



Payne||1982 


| |Narayan et al. 1992; Meszaros 


& Rees||1997; 


Proga et al. |2003 


Uzdensky & MacFadyen 


2006 


1 and the 


magnetic braking 


of the central black hole (Blandford & 


Znajck 1977; Meszaros & Rccs 1997||Lee et al.||2000; McK- 


inney 2006 


Barkov & Komissarov|2008a b 


1. 



A number of groups have studied the potential of mag- 
netic mechanism in the collapsar scenario using Newtonian 
MHD codes and implementing the Paczynski-Witta poten- 
tial in order to approximate the gravitational field of cen- 
tral BH (|Proga et al.|[2003] |Fujimoto et~aT^2006"| |Nagataki 



et al. 



2007) 



In this approach it is impossible to capture 
the Blandford-Znajek effect and only the magnetic braking 
of accretion disk can be investigated. The general conclu- 
sion of these studies is that the accretion disk can launch 
magnetically-driven jets provided the magnetic field in the 
progenitor core is sufficiently strong. This field is further am- 
plified in the disk, partly due to simple winding of poloidal 
component and partly due to the magneto-rotational insta- 
bility (MRI), until the magnetic pressure becomes very large 
and pushes out the surface layers of the disk. In particular, 



Proga et al. (20031 studied the collapse of a star with a 
similar structure to that considered MacFadyen fc Woosley| 
( 1999 1 and with purely monopole initial magnetic field of 
strength B — 2 x 10 14 G at r = 3r 9 . The corresponding to- 
tal magnetic flux, ~ 10 26 Gcm 2 , is in fact within the 
observational constraints on the field of magnetic stars, in- 
cluding white dwarfs ( Schmidt et al.| |2003), Ap-stars (Fer- 



rario & Wickramasinghc 2005J and massive O-stars [Donati 



et al. (2002|. They used realistic equation of state and in- 
cluded the neutrino cooling but not the heating. The simu- 
lations lasted up to t ~ 0.28s in physical time during which 
a Poynting-dominated polar jet has developed in the solu- 
tion. The jet power seemed to show strong systematic de- 
cline, from ~ 10 51 ergs _1 down to ~ 10 50 ergs _1 or even less, 
which might signal transient phenomenon. 



Sekiguchi & Shibata (20071 studied the collapse of ro- 



1 However, they have simulated only the very initial stages of 
the collapsar evolution, t < 2.2s, whereas the jet simulation of 
|MacFadyen fc W oosley ( 1999J started at t = 7s when the density 
in the polar regions above the black hole has dropped down to 
10 6 gcm -3 leading to much lower cooling rates. 



fating stellar cores and formation of BHs and their disks in 
the collapsar scenario using the full GRMHD approxima- 
tion. Their simulations show powerful explosions soon after 
the accretion disk is formed and the free falling plasma of 
collapsing star begins to collide with this disk. However, 
they did not account for the neutrino cooling and the en- 
ergy losses due to photo-dissociation of atomic nuclei so 
these explosions could be similar in nature to the "suc- 
cessful" prompt explosions of early supernova simulations 
( Bethe||1990| ) . |Mizuno eTaL] p004a|bl ) carried out GRMHD 
simulations in the time-independent space-time of a cen- 
tral BH. Their computational domain did not capture the 
BH ergosphere and thus they could not study the role of 
the Blandford-Znajek effect (|Komissarov||2004al 120081). In 



addition, the energy gains/losses due to the neutrino heat- 
ing/cooling were not included and the equation of state 
(EOS) was a simple polytrope. The simulations run for a 
rather short time, ~ 280r 9 /c where r g = GM/c 2 , and jets 
were formed almost immediately, presumably due to the ex- 
tremely strong initial magnetic field. 

The original theory of the Blandford-Znajek mechanism 
was developed within the framework of magnetodynamics 
(the singular limit of RMHD approximation where the in- 
ertia of plasma particles is fully ignored). The pioneering 
numerical studies of this mechanism have shown that it can 
also operate within the full RMHD approximation, at least 
in the regime where the energy density of matter is much 
smaller compared to that of the electromagnetic field, and 
that it can be captured by modern numerical techniques 



(Komissarov 2004b Koide 2004 McKinney & Gammie 2004 



Nagataki|2009| ). However, no systematic numerical study has 
been carried out in order to establish the exact conditions 
for activation of the BZ mechanism. 

The first simulations demonstrating the potential of the 
BZ mechanism in the collapsar scenario were carried out by 
|McKin ncy ( 2006]) who also used simple polytropic EOS and 
considered purely adiabatic case, but captured the effects of 
black hole ergosphere. The initial configuration included a 
torus with poloidal magnetic and a free-falling stellar enve- 
lope. The results show ultrarelativistic jets with the Lorentz 
factor up to 10 emerging from the black hole indicating that 
the Blandford-Znajek effect may play a key role in the pro- 
duction of GRB jets. Barkov & Komissarov ( 2008a|b I added 
a realistic EOS and included the energy losses due to neu- 
trino emission (assuming optically thin regime) and photo- 
dissociation of nuclei. They also reported the development 
of relativistic outflows powered by the central black hole via 
the Blandford-Znajek mechanism (the magnetic braking of 
accretion disk did not seem to contribute much to the jet 
power though). 

The black hole-driven magnetic stellar explosions and 
relativistic jets observed in these exploratory computer sim- 
ulations invite a deeper analysis of this model. Obviously, 
such an outcome cannot be a generic feature of core-collapse 
supernovae. Indeed, only a very small fraction of SNe Ib/Ic 



seem to produce GRBs (Piran 2005 Woosley & Bloom 



2006 1 . So what are the relevant conditions and can they 
be reproduced as the results of stellar evolution? These are 
the main issues we address in this paper. 
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BLANDFORD-ZNAJEK MECHANISM AND 
GRBS 



The rotational energy of Kerr black hole is 
E tot =M h c 2 f 1 (a) ~ 1.8xl0 54 /i(a) | <■«■;<- 



(1) 



where 



A (a) = 1 



i + VT 



a 2 + a' 



1/2 



Mh is the BH mass and a £ [0, 1) is its dimensionless ro- 
tation parameter. For M bil = 2Mq and a = 0.9 this gives 



the enormous value of E„ 



5 x 10 erg which is about 



fifty times higher than the rotational energy of a millisec- 
ond pulsar and well above of what is needed to explain the 
energy of GRBs and associated hypernovae. Even for a rel- 
atively slowly rotating black hole with a = 0.1 this is still 
a respectable number, 2.3 x 10 51 erg. Moreover, because in 
all versions of the collapsar model the black hole spin is 
aligned with the spin of accretion disk this energy reserver 
is continuously replenished via accretion. Thus, as far as the 
availability of energy is concerned the black hole model looks 
very promising indeed. 

The energy release rate is usually estimated using the 
Blandford-Znajek power for the case of force-free monopole 
magnetosphere 



E BZ — 



-(- 

6c V 



(2) 



where Qh is the angular velocity of the BH and is the 
magnetic flux threading one hemisphere of the BH horizon 
(Here it is assumed that the angular velocity of the BH mag- 
netosphere Q — 0.5Qh-)- This formula is quite accurate not 
only for slowly rotating black holes considered in Blandford 
& Znajek ( 1977 I but also for rapidly rotating BHs ( |Komis-| 
sarov|2001 1 . In application to the collapsar problem it gives 
us the following estimate 



E BZ = 1.4 x 10 51 / 2 (a)*L 7 



M h 

M Q 



ergs 



(3) 



where 



f 2 (a)=a 2 (l + i/lT, 



and *Ph j2 7 = Vth/10 27 Gcm 2 . One can see that the power 
of BZ-mechanism is rather sensitive to black hole's mass 
and magnetic flux. Since, the black hole mass is likely to 
be > 3Mq the observed energetics of hypernovae and long 
GRBs requires ^h,27 — 1. This value is comparable with 
the maximal surface flux observed in magnetic stars, Ap 



stars, magnetic white dwarfs, and magnetars (e.g. Ferrario 
|fc Wickramasinghe|20"05] ). Thus, the magnetic field of GRB 
central engine may well be the original field of the progenitor 
star. 

The estimates i ] l|o[ ) show that braking of black holes 
alone can explain the energetics of GRBs and this is why this 
mechanism is often mentioned in the literature on GRBs. 
However, there is another issue to take into consideration. 
The equation |2| is obtained in the limit where the inertia 
of magnetospheric plasma and, to large degree, its gravita- 
tional attraction towards the black hole are ignored. In con- 
trast, the mass density of plasma in the collapsing star may 



be rather high and has to be taken into account. For exam- 
ple, the magnetohydrodynamic waves may become trapped 
in the accretion flow and unable to propagate outwards. In 
such a case, it would not be possible to extract magnetically 
any of the black hole rotational energy and drive stellar ex- 
plosion irrespectively of how high this energy is. Since the 
BZ mechanism can only operate within the black hole ergo- 
sphere (e.g. Komissarov|2 004a, 2008f, the magnetohydrody- 
namic waves must at least be able to cross the ergosphere in 
the outward direction for its activation. This agrees with the 



conclusion reached in Takahashi et al. ( 1990 1 that the Alfven 



surface of a steady-state ingoing wind must be located in- 
side the ergosphere to have outgoing direction of the total 
energy flux. Thus we should consider the following condition 
for activation of the BZ-mechanism in accretion flows: the 
Alfven speed has to exceed the local free fall speed at the 
ergosphere, 

v a > Vf. 

In order to simplify calculations let us apply the Newtonian 
expressions for both speeds. The Alfven speed is 

« 2 = — 
a 4np ' 

where B is the magnetic field strength and p is the plasma 
mass density, whereas the free-fall speed is 

2GM h 



At r 



2 

V f = 



2r g = 2GM/c the local criticality condition reads 



A-Kpc 



< 1. 



(4) 



That is the energy density of magnetic field has to exceed 
that of plasma in the vicinity the black hole. 

For spherically symmetric flows the condition v a > Vf 
can be written as a constraint on the mass accretion rate, 
M, and the flux of radial magnetic field, ^, 



2nry / Mv f (r) 



> 1. 



Since Vf oc r this condition is bound not to be satisfied 
at large r but we only need to apply it at the ergosphere. 



Using r ~ 2r g we then obtain 



where 



(5) 



(6) 



and k c = 1. In fact, the critical value of k must depend 
on the black hole spin and tend to infinity for small a. The 
Newtonian analysis cannot capture this effect. The relativis- 



tic one, on the other hand, appears rather complicated ( Ca- 
|menzi nd 1989; Takahas hi et al.] |T990). Fortunately, nowa- 
days this issue can be investigated via numerical simula- 
tions. Another, interesting and important issue for time- 
dependent simulations is whether there is only one bifurca- 
tion separating solutions with switched-on and switched-off 
BZ mechanism or the transition is reacher and allows more 
complicated time-dependent solutions, e.g. quasi-periodic or 
chaotic ones. 
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3 TEST SIMULATIONS. 

The first type of simulations presented in this paper are 
designed to test the validity of our arguments behind the 
criterion for activation of BZ-mechanism which we derived 
in Sec[2] using a combination of Newtonian and relativistic 
physics. For this purpose we consider a more or less spheri- 
cal accretion of cold magnetized plasma with vanishing an- 
gular momentum onto a rotating black hole. In order to 
focus on the MHD aspects of the problem we ignore the mi- 
crophysics important in the collapsar problem and consider 
plasma with simplified polytropic EOS. 

The main details of our numerical method and various 



test simulations are described in Komissarov ( 1999 2004b 
[2006] ). The only really new feature here is the introduction 
of HLL-solver which is activated when our linear Riemann 
solver fails. This usually occurs in regions of relativistically 
high magnetization where the magnetic energy density ex- 
ceeds that of matter by more than one order of magnitude 
or in strong rarefactions. This makes the scheme a little bit 
more robust though we are still forced to use a "density 
floor" , that is a lower limit on the value of mass-energy den- 
sity of matteij^] 

Both in the numerical scheme and throughout this sec- 
tion of the paper we utilize geometric units where G = c = 
Mh = 1. The Kerr spacetime is described using the Kerr- 
Schild coordinates so the metric form reads 



ds 2 



where 



gttdt 2 + 2g t 4,dtd<j> + 2g tr dtdr+ 
g^d(f> 2 + 2g r4 ,d(j)dr + g rr dr 2 + geed9 2 , 



(7) 



gu = C- i, 



gu 



-(a sin 6, g tr = <> 



g^ = £sin 2 0/A, g r ^ = — asin 2 0(1 + (), 



and 
A 

C 
£ 
A 



g rr = 1 + C, 



r 2 + a 2 cos 2 d, 
2r/A, 

(r 2 + a 2 ) 2 — a 2 A sin 2 6 
r 2 + a 2 — 2r. 



= A, 



3.1 Setup 

The initial solution has to describe both the flow and the 
magnetic field. In this case the flow is given by the steady- 
state solution for unmagnetised cold plasma (dust). At in- 
finity it is spherically symmetric and has vanishing specific 
angular momentum and radial velocity. The components of 
4- velocity of this flow in the coordinate basis of Kerr-Schild 
coordinates, {d v }, are 



2 We have found that a total switch to the HLL-solver noticeably 
degrades numerical solutions via increasing numerical diffusion 
and dissipation. This is particularly noticeable for slowly evolving 
flow components, like accretion disks (see alsojMignonc fc Bodo| 
1 20061). 



«' = 1 + 07/(1 + 77); 
u* = -a/{A(l + V )); 
u r = -O7; 
u = 0, 



(8) 



where 77 = U (r 2 + a 2 )/2r and its rest mass density is 



p = p+ 



where p+ is the rest mass density at the outer event horizon 
radius (see Appendix |A|) . 

The initial magnetic field is radial with monopole topol- 
ogy 

B 4, = aB r {A + 2r)/(AA + 2r(r 2 + a 2 )); 

B r = Bo sin 0/^7; (10) 

B e = 0. 

The azimuthal component is introduced in order to ensure 
that the electromagnetic fluxes of energy and angular mo- 
mentum vanish (see Appendix[5| and, thus, the initial mag- 
netic configuration is more or less dynamically passive. 

The two-dimensional axisymmetric computational do- 
main is (ro < r < ri) x (0 < 8 < 7r), where tq = 1 + 0.5\/l — fl2 
and ri = 500. Notice that the inner boundary is inside of 
the outer event horizon, r+ = 1 + \/l — a?. This can be 
done because of the non-singular nature of the horizon in 
the Kerr-Schild coordinates. The computational grid is uni- 
form in 0, where it has 100 cells, and almost uniform in log r, 
where it has 170 cells. The metric cell sizes are the same in 
both directions. The grid is divided into rings such that for 
any ring the metric size of its outer cells is twice that of its 
inner cells. The solution on each ring is advanced with its 
own time step, the outer ring time step being twice that of 
the inner ring. At the interface separating cells of different 
rings the conservation of mass, energy, and momentum and 
magnetic flux is ensured via simple averaging procedure. 

At the inner boundary, r = ro, we impose the "non- 
reflective" free-flow boundary conditions. The use of these 
conditions is justified by the fact that this boundary lies in- 
side the outer event horizon in the region where all physical 
waves move inwards. At the outer boundary, r = ri we fix 
the flow parameters corresponding to the initial solution. 
In all runs the simulations are terminated well before any 
strong perturbation generated near the centre reaches this 
boundary. At the polar axis, 6 — 0, n, the boundary condi- 
tions are dictated by the axial symmetry. Namely, quantities 
that do not vanish on the axis are reflected without change 
of sign, whereas quantities that vanish on the axis are re- 
flected with change of sign. 

We consider models which differ only by the rotation 
rate of black hole, a, and the mass accretion rate, M. For 
each value of a, we first deal with the model with re ~ 1. De- 
pending on whether the BZ-mechanism is activated in this 
model or not we either increase or reduce p+ (and hence 
M ) by the factor of three and compute another model. We 
proceed in this fashion until the solution becomes qualita- 
tively different. Then we compute one more model based on 
the mean value of p+ for the previous two models. After 
that we select the closest two models with qualitatively dif- 
ferent solutions and show them on the bifurcation diagram 
summarising the results of our study. 
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Figure 1. Free-fall accretion of cold plasma with zero angular momentum onto a Kerr black hole (a = 0.9); the global structure. The 
colour images show log 10 p in dimensionless units and the contours show the magnetic field lines ( To be more specific the contours show 
the levels of poloidal magnetic flux function, >F At any point (r, 8) this function gives the magnetic flux through the axisymmetric loop 
or radius r passing through this point.). Top panel: model with the activation parameter k = 1.2 at time t = 2000r 3 /e. Middle panel: 
model with k = 1.6 and monopole magnetic field at time t = 1000r 9 /c. Bottom panel: model with k = 1.6 and split-monopole magnetic 
field at time t = 1600r s /c. The unit of length is r g . 
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Figure 2. Free-fall accretion of cold plasma with zero angular momentum onto a Kerr black hole (a = 0.9); the central regions. The top 
panels show the model with k = 1.2 at t = 2000r g /c; the middle panels show the model with k = 1.6 and monopole magnetic field at 
t = 1000rg/c, the bottom panels show the model with k = 1.6 and split-monopole magnetic field at t = 1600r 3 /c. The left panels show 
the rest mass density distribution (p for the top panel and log 10 p for others) and the magnetic field lines. The right panels show the 
ratio of the proper rest mass and magnetic energy densities, log 10 /3 P , and the velocity field (the length of velocity arrow is proportional 
to the square root of speed.). 



© 0000 RAS, MNRAS 000, 000-000 



Activation of the BZ mechanism 7 



3.2 Results 

The simulations do confirm the anticipated bifurcation and 
show that for a wide of parameter a the critical value of k is 
indeed close to unity. For example, the model with a — 0.9 
and k — 1.2 is qualitatively different from the model the 
model with a = 0.9 and k = 1.6. The first model quickly 
settles to a steady state accretion solution (see the top pan- 
els of fig[I]and fig(2|. In fact, on the large-scales, this solution 
is almost indistinguishable from the initial solution. On the 
small scales, near the event horizon, the differences become 
more pronounced. In particular, a stationary accretion shock 
appears just outside of the ergosphere. It is responsible for 
the jump in the rest mass density and the kinks of magnetic 
field lines seen in the top left panel of fig[2] Inside the ergo- 
sphere the magnetization parameter \i is larger than unity 
but only just. The magnetic field lines show no rotation and 
thus the BZ-mechanism remains completely switched-off. 

In contrast, the model with a = 0.9 and k — 1.6 ex- 
hibits powerful bipolar outflow which drives a blast wave 
into the accreting flow (see the middle panel of fig[T|. This 
blast waves overturns the accretion both in the polar and in 
the equatorial direction. The black hole develops highly rar- 
efied magnetosphere (see the middle panels of fig[2| which 
rotates with about half of black hole's angular velocity. The 
shocked plasma of accretion flow is prevented from accretion 
onto the black hole - instead it participates in the large scale 
circulations above and below the equatorial plane. The un- 
physical monopole configuration of magnetic field set for this 
model prevents any escape of magnetic flux from the black 
hole. This ensures operation of the BZ-mechanism even if 
the confining accreting envelope is fully dispersed into the 
surrounding space. 

For a more realistic dipolar configuration one would ex- 
pect at least partial escape of magnetic flux once the outflow 
is developed and thus less power in the outflow. However, 
this is unlikely to effect the value of k c . Indeed, consider 
the case with monopole field where the BZ-mechanism re- 
mains switched off. Then, change the polarity of magnetic 
field lines in the northern hemi-sphere. Because the speeds 
of MHD waves and the magnetic stresses remain invariant 
under change of magnetic polarity this will have no effect 
on the flow. For the very same reason, we expect activa- 
tion of the BZ mechanism for exactly the same value of k in 
both magnetic configurations. In order to verify this conclu- 
sion we computed models with the initial magnetic field of 
split-monopole topology and the results fully agree with the 
expectations - the critical value of k remains unaffected but 
the power of BZ-mechanism is reduced due partial escape 
of magnetic flux from the black hole. The bottom panel of 
fig[l] shows the large scale flow developed in the model with 
a — 0.9, k = 1.6, and split-monopole initial magnetic field. 
The flow pattern is very similar to that of the model with 
monopole field but the expansion rate is noticeable lower, 
approximately by the factor of 1.8 in the polar direction 
and even more in the equatorial direction. The bottom pan- 
els of fig[2]show the inner region of this model. One can see 
that the accretion is no longer fully halted - the accreting 
plasma finds its way into the black hole in the equatorial 
zone and the black hole magnetosphere is confined only to 
the polar region of the funnel shape. 

The results for a = 0.9 show that not only there is 
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Figure 3. The bifurcation diagram. The stars show the models 
with the lowest value of K which still show activation of the BZ 
mechanism. The circles show models with the highest value of k 
where the BZ mechanism remains turned off. The pair of points 
with a = 0.9 in the right bottom corner of the diagram correspond 
to proper collapsar simulations. The other points show the results 
for test simulations. 



a bifurcation between the regimes with switched-on and 
switched-off BZ process but also that the critical value of 
the control parameter k, and hence j3 p , is indeed close to 
unity at least for rapidly rotating black holes. The results of 
our parametric study, summarised in fig[3j show that in fact 
K c (a) is a relatively weak function of a. It is confined within 
the range (1, 2) for 0.2 < a < 1.0. 



4 COLLAPSAR SIMULATIONS 

The test simulations described in Sec|3] allow us to deter- 
mine the basic criteria for activation of the BZ mechanism 
in principle. However, their set-up is rather artificial as in 
many astrophysical applications the accreting plasma has 
sufficient angular momentum to ensure formation of accre- 
tion disk and thus to render the approximation of spherical 
accretion unsuitable. This anisotropy of mass inflow is likely 
to make the integral activation condition |5| less helpful. On 
the other hand, one would expect the local condition con- 
dition Q to be more or less robust. In order to check this 
we analysed the results obtained in our current numerical 
study of magnetic collapsar model. 



4.1 Computational grid and boundary conditions 

Like in the test simulations we use the Kerr metric in 
Kerr-Schild coordinates, {(f>, r, 0}. We selected M h = 3Mq 



© 0000 RAS, MNRAS 000, 000-000 



8 S.S. Komissarov & M.V.Barkov 



and the rather optimistic value of a = 0.9 for the cen- 
tral black hole. The two-dimensional axisymmetric compu- 
tational domain is (ro < r < ri) x (0 < 6 < tt), with ro = 
(1 + 0.5V1 ~ a 2 ) r g = 5.4 km and n = 5700 r g = 25000 km. 
The computational grid is of the same type as in the test 
simulations (Sec 3.1 1 but with more cells, 450 x 180. The 



boundary conditions are also the same as in the test simu- 
lations. 



4.2 Microphysics 

For these simulations we use realistic equation of state 
(EOS) that takes into account the contributions from radi- 
ation, lepton gas including pair plasma, and non-degenerate 
nuclei (hydrogen, helium, and oxygen). This is achieved 



via incorporation of the EOS code HELM (Timmes & 
Swcsty 2000), which can be downloaded from the web-site 



http://www.cococubed.com/code_pages7eos.shtmll 



The neutrino cooling is computed assuming opti- 
cally thin regime and takes into account URCA-processes 



(Ivanova et al. 1969), pair annihilation, photo-production, 
and plasma emission ( Schinder et al.|1987 l, as well as syn- 
chrotron neutrino emission | |Bezchastnov et al.|[T997[ ). In 
fact, URCA-processes strongly dominate over other mech- 
anisms in this problem. Photo-disintegration of nuclei is 
included via modification of EOS following the prescrip- 



tion given in Ardeljan et al. (20051. The equation for mass 



fraction of free nucleons is adopted from |Woosley fc Baron| 
( |1992[ ). We have not included the radiative heating due to 
annihilation of neutrinos and antineutrinos produced in the 
accretion disk mainly because this requires elaborate and 
time consuming calculations of neutrino transport. 



4.3 Model of collapsing star 

The collapsing star is introduced using the simple free-fall 
model by Bethe ( 19901. In this model it is assumed that im- 



mediately before the collapse the mass distribution in the 
star satisfies the law pr 3 =const and that the infall pro- 
ceeds with the free-fall speed in the gravitational field of the 
central black hole. This gives us the following radial velocity 

(2GM/r) 1/2 



and mass density 



pCiX 10' 



10 7 cm 



-3/2 



g cm 



(11) 



(12) 



for the collapsing star, where Ci is a constant that depends 
on the star mass |Bethe 1990) and t s is the time since the 
start of collapse. The corresponding accretion rate and ram 
pressure are 



M = 0.056 d 



M 



3M, 



© 



1/2 



M. 



s 



Pram = 



7.9x10^^- (tt-Y'fr) 
3M U s) \10 7 cmJ 



-5/2 



(13) 



,(14) 



respectively. For a core of radius r c — 10 9 cm and mass 
3.0A_Q the collapse duration can be estimated as t c — 
2r c /3->fj ~ 0.6s. Since, in this study we explore the pos- 
sibility of early explosion, that is soon after the core col- 



lapse, we set t a = Is. Because GRBs are currently associ- 
ated with more massive progenitors we consider C\ = 3, 9. 
This gives us the accretion rates of M ~ 0.166A/qs _1 and 
M ~ 0.5-Mqs -1 respectively. 

On top of this we endow the free-falling plasma with 
angular momentum and poloidal magnetic field. The angular 
momentum distribution is 



I = sint 



loir sin (9 / n f 
lo 



where n = 6300km and Iq 



if 
if 

10 1 



rsinS < ri 
rsinO > n 



(15) 



-i 



4.4 Initial magnetic field 



The origin of magnetic field is probably the most impor- 
tant and difficult issue of the magnetic model of GRBs. One 
possibility is that this field already exists in the progeni- 
tor and during the collapse it is simply advected onto the 
black hole. The recent numerical studies indicate that the 
poloidal and the azimuthal components of the progenitor 
field should have similar strengths (Braithwaite & Spruit 



2004). The magnetic flux freezing argument then suggests 



that during the collapse the poloidal field becomes dominant 
and thus in our simulations we may ignore the initial az- 
imuthal field. As to the initial poloidal component we adopt 
the solution for a uniformly magnetized sphere in vacuunj^] 
To be more specific, we assume that inside the sphere of 
radius r m = 4500 km the magnetic field is uniform and has 
strength of either B = 3 x 10 9 , 10 10 , or 3 x 10 10 G, whereas 
outside it is described by the solution for magnetic dipole. 
The selected values for the field strength allow to capture 
the BZ bifurcation. Given the limitations of axisymmetry 
we can only consider the case of aligned magnetic and rota- 
tional axes. 



4.5 Results 

Table [l] describes the key parameters of six models inves- 
tigated in this study and also summarises the outcome of 
simulations. One of the parameters, vp™ 1 "^ [ s the total mag- 
netic flux enclosed by the equator of the uniformly magne- 
tized sphere in the initial solution. Obviously, this gives us 
the highest magnetic flux that can be accumulated by the 
black hole in the simulations. Whether this maximum value 
is actually reached at some point depends on the progen- 
itor rotation. For slow rotation the whole of the uniform 
magnetized core can be swallowed by the hole before the 
development of accretion disk. At this point equals to 
\$i max and then it begins to decline as the result of advec- 
tion of the oppositely directed parts of dipolar loops (These 
loops are clearly seen in fig|5]). For fast rotation the accretion 
disk forms before reaches \^ max an( i after this point the 

3 The intention was to set up a solid body rotation within the 
cylindrical radius r;, however, by mistake an additional factor, 
sin 0, was introduced in the start-up models (The same mistake 



was made in the simulations described in IBarkov fc Komissarovl 
(2008a).) Although somewhat embarrassing, this does not affect 
the main results of this study. 

4 Thus, we ignore the magnetic field that is already threading 
the black hole by the start of simulations. 
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Figure 4. Model B close to the explosion, at t = 237ms. Top left: log 10 p and the poloidal velocity; Top right: the angular velocity of 
plasma and the poloidal velocity; Bottom left: the ratio of azimuthal and poloidal components of magnetic field, log l0 ( Tip / ): an d 
the magnetic field lines; Bottom right: log 10 fl p , the magnetic field lines, and the poloidal velocity field. The length of velocity arrows is 
proportional to the square root of speed. 
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Ci *™ 



Exp * h , 27 E 51 



A 


3 


3 


19.1 


1.07 


Yes 


11.6 


9.6 


B 


3 


9 


19.1 


0.62 


Yes 


13.1 


12.7 


C 


1 


3 


6.67 


0.36 


Yes 


2.80 


0.4 


D 


1 


9 


6.67 


0.21 


No 






E 


0.3 


3 


1.91 


0.11 


No 






F 


0.3 


9 


1.91 


0.06 


No 







Table 1. Summary of numerical models in collapsar simulations. 
Bq,w is the initial magnetic field strength in units of 10 10 G, Ci 
is the density parameter in eq |12| ty%! r ax is the maximum value of 
magnetic flux that can be accumulated by the black hole in units 
of 10 27 Gcm 2 , K ma x is the value of n corresponding to 9^ ax 7 
"Exp" is the explosion indicator, is the actual magnetic 

flux accumulated by the black hole by the time of explosion, E 51 
is the mean power of the BZ mechanism in units of 10 51 erg/s. 



growth of slows down noticeably as the disk accretion is 
much slower compared to the free-fall accretion. Once the 
oppositely directed parts of magnetic loops begin to cross 
the accretion shock and enter the turbulent domain that sur- 
rounds the disk the annihilation of poloidal magnetic field in 
this domain becomes an important factor. As the result the 
maximum value of magnetic flux available in the solution 
begins to decrease before ^ reaches ip max , 

Umax is the value of parameter k computed using \fy max 
instead of the actual value of black hole flux in eq[5] Quick 
inspection shows that the critical value of k based on the 
values of K max is about k c = 0.3 (see models C and D in 
table [I]), almost five time lower than that in the test sim- 
ulations of spherical collapse. In fact, the actual value of 
magnetic flux accumulated by the black hole, 11^, by the 
time of explosion is always noticeably lower than n ma x (see 
table [l]). In model C it is 2.4 times lower thus driving the 
actual critical value of k down to the value around 0.1. This 
result is undoubtedly connected to the highly anisotropic 
nature of the flow in the close vicinity of the black hole 
that develops in the simulations. Most importantly, there is 
a clear separation of the mass and magnetic fluxes: while a 
large fraction of mass is accreted through the accretion disk, 
the bulk of magnetic field lines connected to the black hole 
occupy the space above (and below) the disc (see fig|4|. 

All models considered in this study first pass through 
the passive phase of almost spherical accretion. During this 
phase the black hole accumulates magnetic flux but as 
^max j g not su fg c ; en tiy high the Blandford-Znajek mech- 
anism remains switched-off. Then, as the angular momen- 
tum of plasma approaching the black hole exceeds that of 
the marginally bound orbit, the accretion disk forms in the 
equatorial plane and an accretion shock separates from its 
surface. This shock is best seen in the density plot of fig[4] 
where it appears as a sharp circular boundary between the 
light blue and dark blue areas of radius ~ 30r 9 ~ 135km. In 
fact, the accretion shock oscillates in a manner reminding 
the oscillations found earlier in supernova simulations (e.g. 
Blondin et aT1|2003| |Buras et al-1|2006| |Scheck et aiT1|2008| ) , 
as well as in simulations of accretion flows onto black holes 
(Nagakura & Yamada 2009 1, and attributed to the so-called 
Stationary Accretion Shock Instability (SASI). Most likely, 
the physical origin of shock oscillations in our simulations is 
the same but we have not explored this issue in details. 



During this second phase the differential rotation re- 
sults in winding up of the magnetic field in the disk and its 
corona where the azimuthal component of magnetic field be- 
gins to dominate (top-right and bottom-left panels of fig|4|. 
The corona acts as a "magnetic shield" deflecting the plasma 
coming at intermediate polar angles towards either the equa- 
torial plane or the polar axis (see fig[4]|. In models D,E and 
F the second phase continues till the end of run and the 
Blandford-Znajek process remains switch-off permanently. 
There are no clear indications suggesting that this phase 
may end soon although the accretion disk grows monotoni- 
cally in radius and this may eventually bring about a bifur- 
cation. 

In models A,B,C after several oscillations the solution 
enters the third phase which is characterized by a non-stop 
expansion of the accretion shock which then turns into a 
blast wave. During this transition the black hole develops 
a low mass density magnetosphere which occupies the ax- 
ial funnel through which the rotational energy of the black 
hole is transported away away in the form of Poynting flux 
and powers the blast wave. The typical large-scale structure 
of the the solution during the third phase is shown in fig[5] 
(see also Barkov & Komissarov (2008a) ). It is rather simi- 
lar to the structure observed in the simulations of spherical 
accretion in the case of split-monopole magnetic field (see 
Sec[3~2l). 



Figure [4] shows the solution for model B at time t — 
237ms, which is close to the end of the second phase, and 
helps to understand how exactly the BZ mechanism is ac- 
tivated in our simulations. At this stage, the accretion disk 
is well developed and accounts for ~ 50 — 60% of mass flux 
through the event horizon whereas most of the magnetic 
field lines threading the horizon avoid both the disk and 
corona that surrounds the disk and contains predominantly 
azimuthal magnetic field. The strong turbulent motion in 
the disk and its corona must be the reason for this expul- 
sion of mag netic flux ( |Zerdovich|1957| |Radler||1968| |Tao et 



[aT]p98| . In any case, the observed separation of magnetic 
and mass fluxes results in a rather inhomogeneous magne- 
tization of plasma near the black hole. Most importantly, 
there are regions that have much higher local magnetization 
compared to the case of spherical accretion with the same 
value of integral control parameter k. 

Moreover, the magnetic shield action of the disk corona 
can promote further local increase of magnetization. Indeed, 
the corona can deflect plasma passing through the accretion 
shock at intermediate latitudes in the direction which does 
not necessarily coincide with the local direction of the mag- 
netic field. For example, in figure [4] the stagnation point 
of post-shock flow in the upper hemisphere is located at 
(r, z) ~ (25, 20)r 9 whereas the separation between the mag- 
netic field lines connected to the hole and to the disk occurs 
at much lower latitude, around (r,z) — (28, 10)r 9 . Thus, 
the direct mass loading of the field lines entering the shock 
between these points with is currently suspended. On the 
other hand, at the other end of these field lines, where they 
penetrate the event horizon, their unloading continues as 
usual. As the result, the magnetization along these lines 
increases allowing to reach the local criticality condition, 
f3 p < 1 (see the intense blue regions in the left panel of fig|4|. 
Once this condition is satisfied the BZ mechanism turns on, 
at this point only very locally and the initial rotation rate 
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Figure 5. Model B soon after the explosion. The figure shows 
the solution at t = 332ms. The colour image gives log 10 p and the 
contours show the magnetic field lines. 



of the field lines is significantly lower then 0.5f2h. However, 
an additional magnetic energy is now being pumped into 
the shield via the BZ mechanism promoting its expansion 
and increasing its efficiency as a flow deflector. This causes 
further increase of the magnetization at its base, further in- 
crease of magnetosphere's rotation rate and thus higher BZ 
power and so on. The process develops in a runaway fashion. 

Figure [6] shows the magnetization parameter /3 P at r — 
2r g as a function of the polar angle for model B (at the same 
time as in figQ and for model D at a similar time. One can 
see that for model D the magnetization is lower. In fact, the 
local criticality condition, (5 P < 1, is not satisfied anywhere 




Figure 6. The magnetization parameter f3 p = 4irpc 2 /B 2 at r = 
2r g for models B at t = 237ms (solid line) and D at t = 360ms 
(dashed line). 



for this model and this explains why the explosion is not 
produced in this modej^] 



5 DISCUSSION 



5.1 



The strength and origin of magnetic field. 



There seems to be only two possible origins for the magnetic 
field in the collapsar problem. Either this is a fossil field of 
the progenitor star of the field generated in the accretion 
disk and its corona. Here we consider both possibilities. 

The blunt application of the integral condition (|5| for 
activation of the Blandford-Znajek mechanism gives us the 
critical value of magnetic flux accumulated by the black hole 



4.3 x 1(T 



3M<7> 



Ms 



1/2 



G cm 



Although the strength of magnetic field in the interior of 
massive stars which end up as failed supernovae is not known 
this is uncomfortably high compared to fluxes deduced from 
the observations of "magnetic stars" and rises doubts about 
the fossil hypothesis. However, the actual values of the crit- 
ical magnetic flux found in our simulations are significantly 
lower. As we have already commented in Sec|4] this is due 
to the anisotropy of accretion in the collapsar model. This 
anisotropy is expected to increase with time as more distant 
parts of the progenitor with higher angular momentum enter 



5 These results also suggest that both numerical diffusion and 
resistivity reduce chances of successful explosion in computer sim- 
ulations via smoothing the distributions of mass and magnetic 
field and thus decreasing inhomogeneity of j3 p . 
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the accretion shock (we could not study this late phase of the 
collapse, t > Is, because of various technical issues.). Due 
to the higher angular momentum the predominant motion 
downstream of the shock will be towards the equator and 
this should result in (i) further reduction of the mass flux 
and the ram pressure of accreting plasma in the polar region 
and (ii) more effective mass-unloading of magnetic field lines 



entering the black hole (see Sec 4.5 1. Both effects promote 
higher plasma magnetization in the polar region around the 
black hole which may eventually exceed the local criticality 
condition (jZJ even in models with < 10 27 Gcm 2 , making 
the fossil hypothesis more plausible. 

Our results seem at odds with the Newtonian simula- 
tions by |Proga~e t al. (20031 who obtain magnetically driven 
explosions powered by the disk wind already within 300 ms 
in physical time since the start of simulations for a model 
with similar mass accretion rate, M ~ 0.29Mq/s, but much 
lower magnetic field, \& ~ 10 26 Gcm 2 . Since their initial 



model is similar to that of MacFadyen & Wooslcy ( 1999 1 the 



difference cannot be attributed to the development of highly 
rarefied polar channel seen in the simulations of Mac Fadyen| 
& Woosley ( 1999[ ) as this occurs only several seconds later. 
Other possible reasons include differences in the framework, 
details of the initial solutions, resolution, implementation of 
microphysic^] etc. Hopefully, simulations by other groups 
will soon clarify this issue. 

Another factor that is likely to reduce the constraint 
on the strength of the fossil field, and the one that is com- 
pletely ignored in our simulations, is the plasma heating due 
to neutrino-antineutrino annihilation. In the most extreme 
scenario this heating completely overpowers accretion in the 
polar direction and drives a relativistic GRB jet ( |MacFadyen| 
& Woosley 19991. Although this is an attractive model of 



GRB explosions by itself we also note that the creation of 
low density funnel in this model helps to activate the BZ 
mechanism along the magnetic field lines that are advected 
inside this funnel. Thus, the neutrino heating and BZ mech- 
anism may work hand in hand in driving GRBs and hyper- 
novae. 

Even in failed supernovae the stellar core collapse does 
not lead directly to formation of a black hole but still pro- 
ceeds through the proto-neutron star (PNS) phase first. If 
the period of PNS is very short, around several millisecons, 
its rotational energy is sufficient to strongly modify the di- 
rection of collapse and the initial conditions for the black 
hole phase. In particular, the magnetic stresses due to ei- 
ther fossil or MRI-generated magnetic field can drive pow- 
erful bipolar outflows ( Moisee nko et al.|2006 Obergaulingcr 
et al.||2006| |Burrows et al. | |2007 i. As reported by Burrows 



et al. | ([2007j) these outflows can coexist with equatorial ac- 
cretion and may not be able to prevent eventual collapse of 
PNS. However, they certainly overcome accretion in the po- 
lar direction and create favorable conditions for activation 
of the BZ mechanism. The very fast rotation of progenitor's 
Fe core required in this scenario, with the period around few 
seconds, is unlikely to be achieved in stars with strong fossil 



magnetic field (see Sec 5.2 1 



6 In particular, |Proga et al.| | |2003| state that neutrino cooling 
is dominated by pair annihilation whereas we find that it is the 
lepton capture on free nucleons which is dominant in our models. 



A related issue is whether a significant fraction of fos- 
sil magnetic flux available in the progenitor can actually be 
accumulated by the black hole. The traditional viewpoint 
is that the inward advection of magnetic field in accretion 
disk is effective and can result in build-up of poloidal mag- 
netic field up to the equipartition with the disk gas pres- 



sure ( Bisnovatyi-Kogan & Ruzmaikin 


1974 




1976 


Bland- 


ford & Znajek 1977| |Macdonald & Thorne|11982; Bcgelman 


et al.||1984l. On the other hand, it has also been 


argued 



that in standard viscosity driven accretion disk the turbu- 
lent outward diffusion of poloidal field may easily neutralise 
the inward advection well before it reaches the equipartition 
strength(|Lubow et al.|1994||Heyvaerts et al.|1996||Ghosh fe 



Abramowicz| 1 99*7 Livio et al.|1999 |. However, the disk an 
gular momentum can also be removed by magnetic stresses 
via interaction with the disk corona and wind and this may 
change the balance in favour of the inward accretion of 



poloidal magnetic flux ( Spruit & Uzdensky 20051 Robstein 



& Lovelace 20081. This has been demonstrated in recent 
3D Newtonian simulations of radiatively inefficient accretion 
disks where the poloidal magnetic field reached the equipar- 



tition strength near the central "black hole" ( Igumenshchev 
2008}. In our simulations, a significant fraction of black 
hole's magnetic flux is accumulated before the development 
of accretion disk. However, the flux does not escape after 
the disk is formed. On the contrary, it keeps increasing until 
the oppositely directed magnetic field lines begin to enter 
the accretion shock. Thus, our results are in agreement with 



those by Igumenshchev (2008) and suggest strong magnetic 



interaction between the disk and its corona/wind. 

The alternative origin of magnetic field is via dynamo 
process in the turbulent accretion disk. In this case the sat- 
urated equilibrium magnetic field strength can be estimated 
from the energy equipartition argument 



B? 



(16) 



where Bt and vt are the turbulent magnetic field and velocity 
respectively. However, the typical scale of this field is rather 
small, of the order of the disk hight, h <C r s . When magnetic 
loops of such small scale are advected onto the black hole one 
would expect quick, on the time scale of r g /c, reconnection 
and annihilation of magnetic field, resulting in low residual 
magnetic flux penetrating the black hole. |Tout fe Pringlej 
(1996;) argued that similar reconnection processes in the disk 
corona may lead to inverse cascade of magnetic energy on 
scales exceeding h. The expected spectrum in this range is 
B\ oc 1/A and thus on the scale of r g one could find the 
regular magnetic field of strength 



B„ 



-B t 



The standard theory of a-disks gives 
v t ~ afifc/i, 

where is the Keplerian frequency, h is the half-thickness of 
the disk and a = 0. l-i-0.01 is a free parameter of this theory. 
According to the theory of thin neutrino cooled accretion 



disks by Popham et al. ( 1999 1 



0.09 



Mh 
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and 



-1.7/ \-2.55/ , \ 

ift \ / r \ / Ms \ g 



M / cm 3 ' 



Putting all these equations together we find the typical black 
hole magnetic flux 



* h = 3 x 10" 



26 f a \ 55 f M h 

Mr, 



Ms \ 2 
— — Gem . 

Mo 



The corresponding power of BZ mechanism is 



1.3 x 10 50 /2(a) 



Ms 



erg s 



(17) 



(18) 



almost independently on the black hole mass. One can see 
that even for maximally rotating black holes one should not 
expect much more than 10 51 erg from a collapse of a IOMq 
progenitor. On the other hand, in this model one should ex- 
pect a much more powerful magnetically driven wind from 
the accretion disk ( |Livio et al.||1999| ). Mass loading of this 
wind is likely to be much higher due to the processes similar 
to solar coronal mass injections are the terminal Lorentz fac- 
tors are not expected to be as high as required from the GRB 
observations. However, the disk wind could be responsible 
for the ~ 10 52 ergs required to drive hypernova. Having said 
that, we need to warn that the above equations are purely 
Newtonian and thus their application to the very vicinity of 
black hole may introduce large errors. Moreover, the current 
theory of magnetic dynamo in accretion disks can hardly be 
called well developed and so one could expect few surprises 
as our understanding of the dynamo improves. Finally, the 
above picture is in conflict with our numerical models, par- 
ticularly those which do not show BZ driven explosions. It 
is difficult to pin down the exact reasons for this but the 
key issues are likely to be the condition of axisymmetry and 
numerical resolution. 



5.2 The nature of rotation 

With a typical rotation rate of 200km/s on the zero-age 
main-sequence (e.g. Penny 1996 Howarth et al. 19971 the 
standard evolutionary models of solitary stars predict very 
rapid rotation of stellar cores by the time of collapse ( |Heger| 
et al.||2000| |Hirschi et aLl|2004| |2005[). According to these 



models the typical periods of newly born pulsars in suc- 
cessful supernovae are around one millisecond and "failed 
supernovae" conceive rapidly rotating black holes with ac- 
cretion disks, as required for GRBs. These results, however, 
contradict to observations which clearly indicate that young 
millisecond pulsars and GRBs are much more rare. 

When the magnetic torque due to the relatively weak 
magnetic field generated via the Tayler-Spruit dynamo 



( Tayler 1973 Spruit 2002 I is included in the evolutionary 



models they result in significant spin down of stellar cores 
during the red giant phase and during the intensive mass 
loss period characteristic for massive stars at the Wolf-Rayet 



phase ( Heger et al.|2005 ). As the result the predicted rota- 
tion rates for newly born pulsars come to agreement with 
the observations but on the other hand, the core rotation 
becomes too low compared to what is required for the col- 
lapsar model of GRBs ( [Heger et al.||2005| . Recently, it was 
proposed that low metalicity could be the factor that allows 



to reduce the efficiency of magnetic braking in GRB pro- 
genitors ( |Woosley fc Heger|2006"l |Yoon et al.|2006[ ). On one 
hand, the mass-loss rate of WR-stars decreases significantly 
with metalicity. On the other hand, the evolutionary models 



of Woosley & Heger ( 2006 1 show that low metalicity stars ro- 



tating close to break-up speed may become fully convective, 
chemically homogeneous, and by-pass the red giant phase. 
This result has been confirmed by Yoon et al.| (|2008[) but 



their models also show that even stars with metalicity two 
thousand times below that of Sun inflate above 50 solar radii 
by the stage of carbon burning. For a potential progenitor 
of even long duration GRBs the corresponding light-crossing 
time is uncomfortably high, t lB > 120s. Indeed, this makes 
very difficult to explain the large observed fraction of long 
GRBs with the gamma ray burst duration of only 3-30 sec- 
onds. It is somewhat curious but the high mass loss rates of 
massive stars with solar metalicity has long been considered 
as an important factor in the evolution of GRB progeni- 
tors as it can result in the formation of a compact helium 
star (WR-star) via total loss of extended hydrogen envelop 
I Woosley 1993). Obviously, these complications are not spe- 
cific to the standard fireball model of GRBs. In fact, the 
magnetic field of magnetic stars is much stronger than that 
predicted by the Taylor-Spruit theory of stellar dynamo and 
thus the magnetic braking of GRB progenitors is more se- 
vere in our magnetic model of GRB central engine which 
relies on strong fossil magnetic fieldQ 

To a large extent the above complications are specific 
to solitary stars and can be avoided in close binaries. Firstly, 
the extended envelope of a giant star is only weekly bounded 
by gravity and is easily dispersed in close binary, predomi- 
nantly in the orbital plane ( |Taam fc Sa ndquist 2000). Sec- 
ondly, even in the case of high metalicity and/or strong mag- 
netic field the stellar rotation rate may remain high due to 
the synchronisation of spin and orbital motions. Indeed, the 
minimum separation between the components of binary star, 
L c , is determined by the size of the critical Roche surfaces 
as 



2.MRq" 



(19) 



where R is the radius of more massive component and q 
is the ratio of component masses (q < 1). The rotational 
period is given by 



(1 + 9) 



GM 



(20) 



where L is the actual separation of component and M 
is the mass of more massive component. Using q = 1, 
L — 3R, R — 5R Q and M = 2OM we find that the 
specific angular momentum on the stellar equator is about 
6 x 10 18 cm 2 s _1 . On the other hand the specific angular mo- 
mentum of marginally bound orbit around Schwarzschild 
black hole is 



In 



, GMh 



(21) 



For M h = 2OM this gives us only 3.5 x 10 17 cm 2 s" 1 . Thus, 
when the star collapses one would expect the black hole 



7 The magnetic field of magnetic stars is most likely to be in- 
herited from the interstellar medium during star formation ( |Moss| 



1987| |Braithwaite fc Spruit|2004 1 
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to develop substantial rotation and to get surrounded by 
an accretion disk towards the final phase of the collapse 
I Barkov fc Komissarov|in preparation I. The lower accretion 
rates expected in this phase make the neutrino mechanism 
much less effective but do not have much of effect on the BZ 
mechanism. Indeed, assuming the equipartition of gas and 
magnetic pressure we find that 



8np = 87ra 



-2 2 

pVf 



Applying this at the inner edge (r = r g ) of neutrino-cooled 
accretion disk ( Popham et al. 1999 1 we can estimate the 



maximum magnetic flux that can be kept on the black hole 
by the disk advection as 



3.3 x 10' 



M h 



M Q ) 



Gem 2 . (22) 



For M h = 2OM and M Q = 0.01M Q /s this gives us the 
formidable ~ 7 x 10 2S Gcm 2 . According to eq|3] much 
lower fluxes are needed to explain the observed energetics of 
GRBs. 

Even faster rotation is expected in the case of merger 
of helium star with its compact companion, neutron star or 
a black hole, during the common envelope phase of a close 
binary (jTutukov fc iungelson|1979"||Frayer fc Woosley|1998 



Zhang fc Fryer||2001[ ). Since the orbital angular momentum 
of binary star is 



JUL) 



GMlM 2 L 
M h + M ■ 



(23) 



where M h is the mass of the compact star, we can roughly 
estimate the specific orbital angular momentum deposited 
at the distance r from the centre of the helium star as 



l\r) 



J 2 b (r)/M 2 (r) 



GMlr 



(24) 



M h + M(r) : 

where M(r) is the helium star mass within radius r. When 
M(r) is a slow function, as it is in Bethe's model where 
M(r) 



lnr, 



we may use m eq 



24 



the total mass of the 
helium star, M, in place of M(r). For M bh — 2Mq and 
M hc = 20M Q this gives us 



l(r 9 ) = 1.5 x 10 17 V^cm 2 s 



(25) 

where r g — r/10 9 cm. In spite of all the approximations, such 
estimates agree quite well with the results of numerical sim- 



ulations ( |Zhang fc Fryer 2001 1. For r$ > 1 this is higher than 



the upper value of / = 10 1 ' cm s imposed in the collapsar 



simulations by MacFadyen & Woosley (19991. The higher 



angular momentum is one of the factors that reduces the 
disk accretion rates and hence the efficiency of the neutrino 



annihilation mechanism (Zhang & Fryer 20011 but, as we 



have already commented, for the Blandford-Znajek mecha- 
nism this is much less of a problem. 

The compact companion may initially be a neutron star 
but during the in-spiral it is expected to accumulate substan- 
tial mass and become a black hole ( |Zhang fc Fryer||2001[ ). 
In fact, during the inspiraling the compact companion ac- 
cumulates not only mass but spin as well. The final value 
of the black hole rotation parameter computed in |Zhang fc"| 
Fryer (2001 1 varies from model to model between a — 0.14 



to a — 0.985. Thus, provided that the normal star of the 
common envelope binary is a magnetic star we have all the 
key ingredients for a successful magnetically driven GRB 



explosion: (i) a rapidly rotating black hole, (ii) a rapidly ro- 
tating and (iii) compact collapsing star, and (iv) a "fossil" 
magnetic field which is strong enough to tap the black hole 
energy at rates required by observations (see eq[3]|. It has 
been noted many times that the magnetic fluxes of magnetic 
Ap stars, magnetic white dwarfs, and pulsars are remarkably 



similar (e.g. Ferrario & Wickramasinghe 20051, suggesting 



that they could be related. Now one may consider adding 
another type of objects to this list - GRB collapsars. 



6 CONCLUSIONS 

In this paper we continued our study into the potential of the 
Blandford-Znajek mechanism as main driver in the central 
engines of Gamma Ray Bursts. In particularly, we analysed 
the conditions for activation of the mechanism in the col- 
lapsar model where accumulation of magnetic flux by the 
central black hole is likely to be a byproduct of stellar col- 
lapse. 

It appears that the rotating black hole begins to pump 
energy into the surrounding and overpowers accretion only 
when the rest mass energy density of matter drops below 
the energy density of the electromagnetic field near the er- 
gosphere (see eq|4|. Only under this condition the MHD 
waves generated in the ergosphere can transport energy and 
angular momentum away from the black hole. 

In the case of spherical accretion the above criterion can 
be written as a condition on the total magnetic flux accu- 
mulated by the black hole and the total mass accretion rate 
(see eqj5|. However, in the case of nonspherical accretion, 
characteristic for the collapsar model, this integral condi- 
tion becomes less applicable due to the spacial separation of 
magnetic field and accretion flow inside the accretion shock. 
In fact, our numerical simulations show explosions for sig- 
nificantly lower magnetic fluxes compared to those expected 
from the integral criterion. We have not included the neu- 
trino heating in our numerical models. This is another factor 
that can make activation of the BZ mechanism in the col- 
lapsar setup easier. 

Both the general energetic considerations and the con- 
dition for activation of the Blandford-Znajek mechanism re- 
quire the central black hole to accumulate magnetic flux 
comparable to the highest observed values for magnetic 
stars, ^ ~ 10 27 Gcm 2 . Current evolutionary models of mas- 
sive solitary stars indicate that fossil field of such strength 
should slow down the rotation of helium cores below the 
rate required for formation of accretion disks around black 
holes of failed supernovae. The conflict between strong fossil 
magnetic field and rapid stellar rotation required by the BZ 
mechanism may be resolved in the models where the GRB 
progenitor is a component of close binary. In this case the 
high spin of the progenitor can be sustained at the expense 
of the orbital rotation. Moreover, the strong gravitational 
interaction between components of close binary helps to ex- 
plain the compactness of GRB progenitors deduced from the 
observed durations of bursts. 

The magnetic field can also be generated in the accre- 
tion disk and its corona but in this case the magnetic flux 
through the black hole ergosphere is not sufficiently high 
to explain the power of hypernovae by the BZ mechanism 
alone. On the other hand, in this case the disk is expected to 
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drive MHD wind more powerful than the BZ jet. Although 
such a wind is unlikely to reach the high Lorentz factors 
deduced from the observations of GRB jets (due to high 
mass loading) it could be behind the observed energetics of 
hypernovae. 
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APPENDIX A: DUST ACCRETION IN 
KERR-SCHILD COORDINATES 



The problem of dust accretion has been considered in |Misner] 
et al. ( 1973 I. In the coordinate basis of the Boyer-Lindquist 
coordinates, {d/dx v } the components of 4- velocity of a dust 
particle accreting from rest at infinity and having zero an- 
gular momentum are 



l/ = l + C(r 2 + a 2 )/A; 
U+' = 2</A; 
u r ' = ~CV, 
u 9 ' = 0, 



(Al) 



where r\ — U (r 2 + a 2 )/2r. In order to find the 4-velocity 
components in the coordinate basis of Kerr-Schild coordi- 
nate one can simply use the corresponding transformation 
law ( |Komissarov|2004"a[ ) : 

dt = dt' + (2r/A)dr; 
dcf> = d<t>' + (a/A)dr'; 
dr = dr'; 
d6 = d0'. 



(A2) 



This gives us 

u l = i + cv/(i + v); 



-a/(A(l + n)); 



(A3) 



= 0. 



Notice that in both coordinate systems the particles move 
over a conical surface, 9 =const. Their rotation about the 
symmetry axis is due to the inertial frame dragging effect. 

In order to find the density distribution consider the 
steady-state version of the continuity equation: 

d r (s/^gpu r ) = 0, 

where g — — sin 2 8 A 2 is the determinant of the metric tensor 
in the Kerr-Schild (and Boyer-Lindquist) coordinates. This 
shows that 



-gpu r = C(9), 
where C(8) is some arbitrary function and then that 

C(9) 



p(r,9) = 



r 2 + a 2 



At infinity this function does not depend on the polar angle 
only if C oc sin 6 which implies that p does not depend on 6 
for any r. Denoting as p+ the rest mass density at the outer 
event horizon we may now write 



P = P4 



>r+\ 1 
> r J rj 



(A4) 



APPENDIX B: INITIAL MAGNETIC FIELD 

In this section we employ the 3+1 splitting of electrodynam- 



ics described in Komissarov (2004a) where the electromag- 



netic field is represented by four 3- vectors, D,E,B, and H , 
defined via 



B' 



a F 
a 



D' = aF l 
a 



Hi 



. &ijkF 



jk 



(Bl) 
(B2) 

(B3) 
(B4) 



where *_F M " is the Faraday tensor, F^ v is the Maxwell ten- 
sor, djk is the Levi-Civita tensor of space and a is the lapse 
function. In stationary metric, dtg^n = 0, the evolution of 
these vector fields is described by the Maxwell equations for 
electromagnetic field in matter. The effects of curved space 
time are incorporated via the non-Euclidian spatial metric 
and the constitutive equations 



E = aD + /3xB, 
H — aB — f3xD. 



(B5) 
(B6) 



where a is the lapse function and (3 is the shift vector. In 
terms of these vectors the condition of perfect conductivity, 



E = —v x B. 



(B7) 



where v l = dx l /dt is the usual 3- velocity vector. 

Both in the Boyer-Lindquist and the Kerr-Schild coordi- 
nates the poloidal component of Poynting flux, S 1 = —aTl, 
is given by 



S p — E p x if + E$ x Hp 



(B8) 



whereas the the poloidal component of the angular momen- 
tum flux, V = aT, 



0' 



(Hm)B p , 



(B9) 



where m = d/d<f>. 

Suppose that in the Boyer-Lindquist coordinates the 
magnetic field is purely radial, B' g = B'^ = 0. Since (3 \\ m 
eq |B6| then implies that 

h'^ = b; = o. 

Provided this magnetic field is frozen into the flow described 
by eq|Al| we also find that 



171/ / r->/ 

E = -v^xB r . 
Given these results the equations |B9| and |B8| ensure that 
S p = Lp = 0. 

In order to find the corresponding electromagnetic field 
in the Kerr-Schild coordinates we simply apply the transfor- 



mation law A2|. This way we find that 

Sp = Lp = 

and 

»_n _ B r a(A + 2r) 



B 



0, B" 



AA + 2r(r 2 +a 2 ) 



(B10) 



(Bll) 
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The divergence-free condition requires 
d r (VlB r ) = 

and hence 

where 7 is the determinant of the spatial metric. Provided 
that at infinity (where 7 oc sin 2 ff) B r does not depend on 
the polar angle we then have 



B r 



(B12) 
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